Engineering Logical Inflammation Sensing Circuit for Gut Modulation¶
BioCRNPyler Models and Code for Fig 5C, 6BCD¶
Liana Merk
In [1]:
from bioscrape.types import Model
from bioscrape.simulator import py_simulate_model
from biocrnpyler import *
import numpy as np
import pylab as plt
import tqdm
import pandas as pd
import matplotlib as mpl
import holoviews as hv
import bokeh.plotting
import bokeh.io
from bokeh.models import Range1d
from bokeh.layouts import gridplot
bokeh.io.output_notebook()
hv.extension('bokeh')
from bokeh.models import PrintfTickFormatter
import seaborn as sns
sns.set()
%matplotlib inline
Tetrathionate Two Component System¶
In [2]:
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1,
"k2b": 0.3, "k2u": 10,
"k3b": 0.3, "k3u": 10,
"k4b": 1.6, "k4u": 0.016,
"k5b": 0.0001, "k5u":1,
"k6b": 0.0001, "k6u": 0.1,
"k7b": 1, "k7u": 1,
"k8b": 0.0083, "k8u": 5,
"k9b": 0.3, "k9u": 3,
"k10b": 0.1, "k10u":1,
"k11b": 100, "k11u": 1,
"ktx": 0.09, "ktl": 0.05,
"k_dephos": 50, "delta": .9}
Define all species and reactions¶
Some Basic Species¶
In [3]:
# Ribosome
R = Species("R")
# RNA Polymerase
P = Species("P")
Reaction 1: Polymerase binding to P1d¶
In [4]:
# Consitutive Promoter P1d
P1d = Species("P1d")
# Polymerase bound to P1d
P_P1d = Species("P:P1d")
# P + P1d <-> P1d:P
R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])
Transcription 1 (ttrS and ttrR)¶
In [5]:
# ttrS Transcript
ttrS_T = Species("ttrS_T")
# ttrR Transcript
ttrR_T = Species("ttrR_T")
# P1d:P -> ttrS_T + ttrR_T
Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])
Reaction 2 and translation 1: ttrS¶
In [6]:
# ttrS trasncript bound to Ribosome
ttrS_T_R = Species("ttrS_T:R")
# ttrS
ttrS = Species("ttrS")
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])
# ttrS Translation
Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])
Reaction 3 and translation 2: ttrR¶
In [7]:
# ttrR trasncript bound to Ribosome
ttrR_T_R = Species("ttrR_T:R")
# ttrR
ttrR = Species("ttrR")
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])
# ttrR Translation
Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])
Reaction 4: Tetrathionate binding to ttrS, Causing Phosphorylation¶
In [8]:
# Tetrathionate
tt = Species("tt")
# Phosphorylated ttrS
ttrS_P = Species("ttrS_P")
# ttrS + tt <-> ttrS_P + tt
R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])
Reaction 5: ttrR binding to dephosphorylated ttrS¶
In [9]:
# ttrR bound to ttrS
ttrR_ttrS = Species("ttrR_ttrS")
# ttrR + ttrS <-> ttrR:ttrS
R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])
Reaction 6: ttrR binding to phosphorylated ttrS¶
In [10]:
# ttrR bound to phosphorylated ttrS
ttrR_ttrS_P = Species("ttrR_ttrS_P")
# ttrR + ttrS <-> ttrR:ttrS:P
R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])
Reaction 7: Phosphorylation transfer ttrS to ttrR¶
In [11]:
# Phosphorylated ttrR
ttrR_P = Species("ttrR_P")
# ttrR:ttrS:P <-> ttrR:P + ttrS
R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])
Dephosphorylated¶
In [12]:
# ttrR:P -> ttrR
rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])
Reaction 8: Dimerization of phosphorylated ttrR¶
In [13]:
# ttrR:P_dimer
ttrR_P_dimer = Species("ttrR:P_dimer")
# 2ttrR:P <-> ttrR:P_dimer
R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])
Reaction 9: Phosphorylated ttrR dimer activated promoter pttr0, turns to pttr1¶
In [14]:
# Inactive pttr
pttr0 = Species("pttr0")
# Inactive pttr
pttr1 = Species("pttr1")
# ttrR:P_dimer + pttr0 <-> pttr1
R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])
Reaction 10: Polymerase binding to active promoter¶
In [15]:
# Active pttr to polymerase
pttr1_P = Species("pttr1_P")
# P + pttr1 <-> pttr1:P
Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])
Transcription 2: active pttr produces GFP transcript¶
In [16]:
# Active pttr to polymerase
GFP_T = Species("GFP_T")
# P + pttr1 <-> pttr1:P
R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])
Reaction 11: Ribosome binding to GFP transcript and translation 3¶
In [17]:
# Active pttr to polymerase
GFP_T_R = Species("GFP_T_R")
# GFP
GFP = Species("GFP")
# P + pttr1 <-> pttr1:P
R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])
# ttrS Translation
Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])
Degradation Reactions¶
In [18]:
# RNA deg
rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])
# Protein Deg
rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])
Now, pull all the reactions together¶
In [19]:
reactions = [R1, R2, R3, R4, R5, R6, R7, R8, R9, R10, R11,
Rtx1, Rtx2,
Rtl1, Rtl2, Rtl3,
rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
In [20]:
species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
In [21]:
# Create CRN
CRN_tt = ChemicalReactionNetwork(species = species, reactions = reactions)
print(CRN_tt.pretty_print(show_rates = True))
In [22]:
# An initial condition for each species (uninitialized species default to 0)
x0 = {
"tt":200,
"P1d":5,
"pttr0":5,
"P":25,
"R":100,
}
timepoints = np.linspace(0, 600, 100)
In [23]:
Results = CRN_tt.simulate_with_bioscrape(timepoints,x0)
In [24]:
Results.head()
Out[24]:
Holoviews plotting of Bioscrape results¶
Sort out Plotting¶
In [25]:
def plot_holoviews_ttr_bioscrape(results_melted):
plot_group = []
for row in results_melted.iterrows():
sp = row[1]['species']
if sp in ["P1d", "P:P1d", "pttr0", "pttr1", "pttr1_P"]:
plot_group.append('promoters')
elif sp in ["ttrS_T", "ttrR_T", "ttrR_T:R", "ttrS_T:R", "GFP_T", "GFP_T_R"]:
plot_group.append('mRNA')
elif sp in ["ttrS_P", "ttrR_ttrS", "ttrR_ttrS_P", "ttrR_P", "ttrR:P_dimer"]:
plot_group.append('complexes')
elif sp in ["ttrR", "ttrS", "GFP"]:
plot_group.append('proteins')
elif sp in ["P", "R", "tt"]:
plot_group.append('machinery and signals')
else:
print(sp)
results_melted['plot_group'] = plot_group
plot_list = []
for group in results_melted['plot_group'].unique():
plot_list.append(hv.Curve(
data=results_melted.loc[results_melted['plot_group'] == group],
kdims=['time', 'value'],
vdims = ['species']
).groupby(
[ 'species']
).opts(
tools=['hover'],
# legend_position = 'right',
show_grid=True,
title = f'{group}',
width = 500,
#shared_axes=False
).overlay('species').opts(legend_position = 'right',))
return plot_list
Now, proceed with plotting¶
In [26]:
results_melted = pd.melt(Results, var_name='species', id_vars=['time'])
results_melted.head()
Out[26]:
In [27]:
plot_list = plot_holoviews_ttr_bioscrape(results_melted)
hv.Layout(plot_list).cols(2).opts(shared_axes=False)
Out[27]:
Tuning Function¶
In [28]:
def create_reactions(parameters):
# P + P1d <-> P1d:P
R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])
# P1d:P -> ttrS_T + ttrR_T
Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])
# ttrS Translation
Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])
# ttrR Translation
Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])
# ttrS + tt <-> ttrS_P + tt
R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])
# ttrR + ttrS <-> ttrR:ttrS
R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])
# ttrR + ttrS <-> ttrR:ttrS:P
R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])
# ttrR:ttrS:P <-> ttrR:P + ttrS
R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])
# ttrR:P -> ttrR
rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])
# 2ttrR:P <-> ttrR:P_dimer
R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])
# ttrR:P_dimer + pttr0 <-> pttr1
R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])
# P + pttr1 <-> pttr1:P
Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])
# P + pttr1 <-> pttr1:P
R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])
# P + pttr1 <-> pttr1:P
R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])
# ttrS Translation
Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])
# RNA deg
rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])
# Protein Deg
rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])
reactions = [R1, R4, R5, R6, R7, R8, R9, R10, R11,
Rtx1, Rtx2,
Rtl1, Rtl2, Rtl3,
rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
return reactions
RBS Tuning in Tetrathionate Circuit¶
We will change k2u, and k3u, the reaction rates that define ribosome unbinding from ttrS and ttrR transcripts, respectively. We will save the simulations as csv's to plot in another notebook such that the styles match well.
In [29]:
all_sims = pd.DataFrame()
In [30]:
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1,
"k2b": 3, "k2u": 10,
"k3b": 0.00003, "k3u": 10,
"k4b": 5, "k4u": 0.5,
"k5b": 5, "k5u":6,
"k6b": 5, "k6u": 1,
"k7b": 0.01, "k7u": 1,
"k8b": 80, "k8u": 5,
"k9b": 0.004, "k9u": 3,
"k10b": 50, "k10u":1,
"k11b": 10, "k11u": 1,
"ktx": 0.1, "ktl": 1,
"k_dephos": 7, "delta": .9}
In [31]:
def plot_ttr_RBS_tuning(parameters, x0, k2u, k3u, title, timepoints):
experimental_ttr_concs = np.array([ 0, 0.1, 0.25, 0.5, 0.75, 1, 5, 10])
# Get color palette in order
r = list(bokeh.palettes.viridis(10))
r.reverse()
r = r[2:]
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = k2u)
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = k3u)
reactions = create_reactions(parameters)
reactions.append(R2)
reactions.append(R3)
species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
# Recompile the CRN with the right parameters
CRN_curr = ChemicalReactionNetwork(species = species, reactions = reactions)
p = bokeh.plotting.figure(width = 300, height= 250, x_axis_label='Time (hr)', y_axis_label = '[GFP]', title = title)
i = 0
for tt_conc in experimental_ttr_concs:
x0["tt"] = tt_conc
Results = CRN_curr.simulate_with_bioscrape(timepoints,x0)
all_sims[f'{title}_{tt_conc}'] = Results['GFP']
p.line(timepoints/60, Results['GFP'],color = r[i], legend_label = f'{tt_conc}', line_width = 2)
i+=1
p.legend.location = 'top_left'
if k2u != 20:
p.legend.visible = False
p.y_range = Range1d(-0.1,2.5)
return p
In [32]:
timepoints = np.linspace(0, 1400, 1400)
title_list = ['Very Low (LM15) Simulation', 'Low (LM14) Simulation', 'Medium (LM18) Simulation', 'High (LM19) Simulation']
In [33]:
x0 = {
"tt":200,
# PSC101 has copy number 5
"P1d":5,
"pttr0":5,
# Polymerase: 15 uM assuming ~10000 Polymerase molecules / E. coli with a volume 1 um^3
"P":15,
# Ribosome: 120 uM assuming ~72000 Ribosomes / E. coli with a volume 1 um^3
"R":120,
}
In [34]:
parameters = {"k1b": 10, "k1u": 0.0001,
"k2b": 0.3, "k2u": 10,
"k3b": 0.3, "k3u": 10,
"k4b": 1.6, "k4u": 0.016,
"k5b": 0.0001, "k5u":1,
"k6b": 0.0001, "k6u": 0.1,
"k7b": 1, "k7u": 1,
"k8b": 0.0083, "k8u": 0.5,
"k9b": 0.3, "k9u": 3,
"k10b": 10, "k10u":0.0001,
"k11b": 100, "k11u": 1,
"ktx": 0.1, "ktl":0.33,
"k_dephos": 50, "delta": .9}
In [35]:
# Define higher off rates
k2u = [20, 12, 11, 10]
k3u = [15, 8, 6, 5]
plot_list = []
for i in range(len(k3u)):
plot_list.append(plot_ttr_RBS_tuning(parameters, x0, k2u[i], k3u[i], title_list[i], timepoints))
bokeh.io.show(gridplot([plot_list]))
Tetrathionate Two Component System¶
In [2]:
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1,
"k2b": 0.3, "k2u": 10,
"k3b": 0.3, "k3u": 10,
"k4b": 1.6, "k4u": 0.016,
"k5b": 0.0001, "k5u":1,
"k6b": 0.0001, "k6u": 0.1,
"k7b": 1, "k7u": 1,
"k8b": 0.0083, "k8u": 5,
"k9b": 0.3, "k9u": 3,
"k10b": 0.1, "k10u":1,
"k11b": 100, "k11u": 1,
"ktx": 0.09, "ktl": 0.05,
"k_dephos": 50, "delta": .9}
Define all species and reactions¶
Some Basic Species¶
In [3]:
# Ribosome
R = Species("R")
# RNA Polymerase
P = Species("P")
Reaction 1: Polymerase binding to P1d¶
In [4]:
# Consitutive Promoter P1d
P1d = Species("P1d")
# Polymerase bound to P1d
P_P1d = Species("P:P1d")
# P + P1d <-> P1d:P
R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])
Transcription 1 (ttrS and ttrR)¶
In [5]:
# ttrS Transcript
ttrS_T = Species("ttrS_T")
# ttrR Transcript
ttrR_T = Species("ttrR_T")
# P1d:P -> ttrS_T + ttrR_T
Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])
Reaction 2 and translation 1: ttrS¶
In [6]:
# ttrS trasncript bound to Ribosome
ttrS_T_R = Species("ttrS_T:R")
# ttrS
ttrS = Species("ttrS")
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])
# ttrS Translation
Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])
Reaction 3 and translation 2: ttrR¶
In [7]:
# ttrR trasncript bound to Ribosome
ttrR_T_R = Species("ttrR_T:R")
# ttrR
ttrR = Species("ttrR")
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])
# ttrR Translation
Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])
Reaction 4: Tetrathionate binding to ttrS, Causing Phosphorylation¶
In [8]:
# Tetrathionate
tt = Species("tt")
# Phosphorylated ttrS
ttrS_P = Species("ttrS_P")
# ttrS + tt <-> ttrS_P + tt
R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])
Reaction 5: ttrR binding to dephosphorylated ttrS¶
In [9]:
# ttrR bound to ttrS
ttrR_ttrS = Species("ttrR_ttrS")
# ttrR + ttrS <-> ttrR:ttrS
R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])
Reaction 6: ttrR binding to phosphorylated ttrS¶
In [10]:
# ttrR bound to phosphorylated ttrS
ttrR_ttrS_P = Species("ttrR_ttrS_P")
# ttrR + ttrS <-> ttrR:ttrS:P
R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])
Reaction 7: Phosphorylation transfer ttrS to ttrR¶
In [11]:
# Phosphorylated ttrR
ttrR_P = Species("ttrR_P")
# ttrR:ttrS:P <-> ttrR:P + ttrS
R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])
Dephosphorylated¶
In [12]:
# ttrR:P -> ttrR
rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])
Reaction 8: Dimerization of phosphorylated ttrR¶
In [13]:
# ttrR:P_dimer
ttrR_P_dimer = Species("ttrR:P_dimer")
# 2ttrR:P <-> ttrR:P_dimer
R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])
Reaction 9: Phosphorylated ttrR dimer activated promoter pttr0, turns to pttr1¶
In [14]:
# Inactive pttr
pttr0 = Species("pttr0")
# Inactive pttr
pttr1 = Species("pttr1")
# ttrR:P_dimer + pttr0 <-> pttr1
R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])
Reaction 10: Polymerase binding to active promoter¶
In [15]:
# Active pttr to polymerase
pttr1_P = Species("pttr1_P")
# P + pttr1 <-> pttr1:P
Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])
Transcription 2: active pttr produces GFP transcript¶
In [16]:
# Active pttr to polymerase
GFP_T = Species("GFP_T")
# P + pttr1 <-> pttr1:P
R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])
Reaction 11: Ribosome binding to GFP transcript and translation 3¶
In [17]:
# Active pttr to polymerase
GFP_T_R = Species("GFP_T_R")
# GFP
GFP = Species("GFP")
# P + pttr1 <-> pttr1:P
R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])
# ttrS Translation
Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])
Degradation Reactions¶
In [18]:
# RNA deg
rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])
# Protein Deg
rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])
Now, pull all the reactions together¶
In [19]:
reactions = [R1, R2, R3, R4, R5, R6, R7, R8, R9, R10, R11,
Rtx1, Rtx2,
Rtl1, Rtl2, Rtl3,
rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
In [20]:
species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
In [21]:
# Create CRN
CRN_tt = ChemicalReactionNetwork(species = species, reactions = reactions)
print(CRN_tt.pretty_print(show_rates = True))
In [22]:
# An initial condition for each species (uninitialized species default to 0)
x0 = {
"tt":200,
"P1d":5,
"pttr0":5,
"P":25,
"R":100,
}
timepoints = np.linspace(0, 600, 100)
In [23]:
Results = CRN_tt.simulate_with_bioscrape(timepoints,x0)
In [24]:
Results.head()
Out[24]:
Holoviews plotting of Bioscrape results¶
Sort out Plotting¶
In [25]:
def plot_holoviews_ttr_bioscrape(results_melted):
plot_group = []
for row in results_melted.iterrows():
sp = row[1]['species']
if sp in ["P1d", "P:P1d", "pttr0", "pttr1", "pttr1_P"]:
plot_group.append('promoters')
elif sp in ["ttrS_T", "ttrR_T", "ttrR_T:R", "ttrS_T:R", "GFP_T", "GFP_T_R"]:
plot_group.append('mRNA')
elif sp in ["ttrS_P", "ttrR_ttrS", "ttrR_ttrS_P", "ttrR_P", "ttrR:P_dimer"]:
plot_group.append('complexes')
elif sp in ["ttrR", "ttrS", "GFP"]:
plot_group.append('proteins')
elif sp in ["P", "R", "tt"]:
plot_group.append('machinery and signals')
else:
print(sp)
results_melted['plot_group'] = plot_group
plot_list = []
for group in results_melted['plot_group'].unique():
plot_list.append(hv.Curve(
data=results_melted.loc[results_melted['plot_group'] == group],
kdims=['time', 'value'],
vdims = ['species']
).groupby(
[ 'species']
).opts(
tools=['hover'],
# legend_position = 'right',
show_grid=True,
title = f'{group}',
width = 500,
#shared_axes=False
).overlay('species').opts(legend_position = 'right',))
return plot_list
Now, proceed with plotting¶
In [26]:
results_melted = pd.melt(Results, var_name='species', id_vars=['time'])
results_melted.head()
Out[26]:
In [27]:
plot_list = plot_holoviews_ttr_bioscrape(results_melted)
hv.Layout(plot_list).cols(2).opts(shared_axes=False)
Out[27]:
Tuning Function¶
In [28]:
def create_reactions(parameters):
# P + P1d <-> P1d:P
R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])
# P1d:P -> ttrS_T + ttrR_T
Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])
# ttrS Translation
Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])
# ttrR Translation
Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])
# ttrS + tt <-> ttrS_P + tt
R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])
# ttrR + ttrS <-> ttrR:ttrS
R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])
# ttrR + ttrS <-> ttrR:ttrS:P
R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])
# ttrR:ttrS:P <-> ttrR:P + ttrS
R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])
# ttrR:P -> ttrR
rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])
# 2ttrR:P <-> ttrR:P_dimer
R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])
# ttrR:P_dimer + pttr0 <-> pttr1
R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])
# P + pttr1 <-> pttr1:P
Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])
# P + pttr1 <-> pttr1:P
R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])
# P + pttr1 <-> pttr1:P
R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])
# ttrS Translation
Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])
# RNA deg
rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])
# Protein Deg
rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])
reactions = [R1, R4, R5, R6, R7, R8, R9, R10, R11,
Rtx1, Rtx2,
Rtl1, Rtl2, Rtl3,
rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
return reactions
RBS Tuning in Tetrathionate Circuit¶
We will change k2u, and k3u, the reaction rates that define ribosome unbinding from ttrS and ttrR transcripts, respectively. We will save the simulations as csv's to plot in another notebook such that the styles match well.
In [29]:
all_sims = pd.DataFrame()
In [30]:
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1,
"k2b": 3, "k2u": 10,
"k3b": 0.00003, "k3u": 10,
"k4b": 5, "k4u": 0.5,
"k5b": 5, "k5u":6,
"k6b": 5, "k6u": 1,
"k7b": 0.01, "k7u": 1,
"k8b": 80, "k8u": 5,
"k9b": 0.004, "k9u": 3,
"k10b": 50, "k10u":1,
"k11b": 10, "k11u": 1,
"ktx": 0.1, "ktl": 1,
"k_dephos": 7, "delta": .9}
In [31]:
def plot_ttr_RBS_tuning(parameters, x0, k2u, k3u, title, timepoints):
experimental_ttr_concs = np.array([ 0, 0.1, 0.25, 0.5, 0.75, 1, 5, 10])
# Get color palette in order
r = list(bokeh.palettes.viridis(10))
r.reverse()
r = r[2:]
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = k2u)
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = k3u)
reactions = create_reactions(parameters)
reactions.append(R2)
reactions.append(R3)
species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
# Recompile the CRN with the right parameters
CRN_curr = ChemicalReactionNetwork(species = species, reactions = reactions)
p = bokeh.plotting.figure(width = 300, height= 250, x_axis_label='Time (hr)', y_axis_label = '[GFP]', title = title)
i = 0
for tt_conc in experimental_ttr_concs:
x0["tt"] = tt_conc
Results = CRN_curr.simulate_with_bioscrape(timepoints,x0)
all_sims[f'{title}_{tt_conc}'] = Results['GFP']
p.line(timepoints/60, Results['GFP'],color = r[i], legend_label = f'{tt_conc}', line_width = 2)
i+=1
p.legend.location = 'top_left'
if k2u != 20:
p.legend.visible = False
p.y_range = Range1d(-0.1,2.5)
return p
In [32]:
timepoints = np.linspace(0, 1400, 1400)
title_list = ['Very Low (LM15) Simulation', 'Low (LM14) Simulation', 'Medium (LM18) Simulation', 'High (LM19) Simulation']
In [33]:
x0 = {
"tt":200,
# PSC101 has copy number 5
"P1d":5,
"pttr0":5,
# Polymerase: 15 uM assuming ~10000 Polymerase molecules / E. coli with a volume 1 um^3
"P":15,
# Ribosome: 120 uM assuming ~72000 Ribosomes / E. coli with a volume 1 um^3
"R":120,
}
In [34]:
parameters = {"k1b": 10, "k1u": 0.0001,
"k2b": 0.3, "k2u": 10,
"k3b": 0.3, "k3u": 10,
"k4b": 1.6, "k4u": 0.016,
"k5b": 0.0001, "k5u":1,
"k6b": 0.0001, "k6u": 0.1,
"k7b": 1, "k7u": 1,
"k8b": 0.0083, "k8u": 0.5,
"k9b": 0.3, "k9u": 3,
"k10b": 10, "k10u":0.0001,
"k11b": 100, "k11u": 1,
"ktx": 0.1, "ktl":0.33,
"k_dephos": 50, "delta": .9}
In [35]:
# Define higher off rates
k2u = [20, 12, 11, 10]
k3u = [15, 8, 6, 5]
plot_list = []
for i in range(len(k3u)):
plot_list.append(plot_ttr_RBS_tuning(parameters, x0, k2u[i], k3u[i], title_list[i], timepoints))
bokeh.io.show(gridplot([plot_list]))
In [36]:
all_sims['Time (hr)'] = timepoints
all_sims.to_csv('./data/all_sims.csv', index = False)
HRP AND Circuit¶
This $\sigma 54$-dependent split activator requires two proteins, hrpR and hrpS, to function. For now, let's just spike in the two proteins and see how our and gate responds.
In [37]:
txtl = CRNLab("GFP")
txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = './parameters/BasicExtract.tsv')
#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
Re1.head()
Out[37]:
Vary hrpR and hrpS proteins¶
Now, we will make sure the AND gate (combinatorial promoter) is working by varying protein concentrations.
Get plotting in order¶
Let's create a plotting function. We will plot a heatmap of protein GFP at the last slice of the simulation. Because we have no degradation of protein, the last timepoint will be the max. We will be able to pass in normalize, which divides all values by the maximum of all simulated values, or max_to_1, which will set the scale bar to be 0 to 1. This basically causes plots that might've been
In [38]:
def simulate_heatmap(inducer_1_name,
inducer_1_values,
inducer_2_name,
inducer_2_values,
title,
x0=x0,
normalize=True,
max_to_1 = False,
CRN=CRN_extract_1,
timepoints=timepoints):
GFP_max = pd.DataFrame()
for conc_1 in inducer_1_values:
x0[inducer_1_name] = conc_1
for conc_2 in inducer_2_values:
x0[inducer_2_name] = conc_2 #Change my initial condition dictionary
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
GFP_max = GFP_max.append({inducer_1_name:conc_1,
inducer_2_name:conc_2,
'GFP_max': Re1["protein_GFP"].values[-1]}, ignore_index=True)
data = pd.pivot_table(data = GFP_max, index = inducer_1_name,
columns = inducer_2_name,
values = 'GFP_max')
fig, ax = plt.subplots(figsize=(6, 4.5))
if normalize:
data /= np.max(GFP_max['GFP_max'])
if max_to_1:
heat_map = sns.heatmap(data,
cmap = 'viridis',
linewidths=.1, vmin = 0, vmax = 1,
cbar_kws={'label': 'Normalized GFP'})
else:
heat_map = sns.heatmap(data,
cmap = 'viridis',
linewidths=.2,
cbar_kws={'label': 'Normalized GFP'})
sns.set(font_scale=1.4)
# plt.tick_params(axis='both', labelsize=20)
plt.xticks(fontsize=18, fontstyle='normal', rotation = 90)
plt.yticks(fontsize=18, fontstyle='normal')
plt.xlabel(inducer_1_name, fontsize=18, fontweight='bold')
plt.ylabel(inducer_2_name, fontsize=18, fontweight='bold')
heat_map.set_title(title)
fig.tight_layout()
return fig, heat_map
And Gate with leak in hrpR¶
Let's say hrpR mutated and now it can recruit RNAP better than it did beter. We will change phrpL_hrpR_RNAP ku = 5 instead of 500.
In [39]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak5_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 1x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_1.png", bbox_inches = "tight")
plt.show()
In [40]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak50_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 10x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_2.png", bbox_inches = "tight")
plt.show()
In [41]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 100x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_3.png", bbox_inches = "tight")
plt.show()
Full Compiled AND gate with expression of hrpR, S¶
LacI promoter class¶
In [42]:
class LacIPromoter2(Promoter):
def __init__(self, name = "pLac", IPTG = "IPTG",
lacI = "lacI", leak = True, assembly = None,
transcript = None, length = 0, mechanisms = {},
parameters = {}, **keywords):
Promoter.__init__(self, name=name, transcript=transcript, assembly=assembly, length=length, parameters=parameters, **keywords)
#RegulatedPromoter.__init__(self = self, name = name, regulators = regulators)
# Transcription is already included in promoter class
self.default_mechanisms = {"binding": One_Step_Cooperative_Binding()}
self.leak = leak
self.IPTG = self.set_species(IPTG, material_type = "inducer")
self.lacI = self.set_species(lacI, material_type = "protein")
self.plac_2_lacI = None
self.IPTG_2_lacI = None
def update_species(self):
mech_tx = self.mechanisms['transcription']
mech_b = self.mechanisms['binding']
species = []
self.complexes = []
"""A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
2lacI <--> 2xlacI
2xlacI + pLac <--> 2xlacI:pLac
"""
species_b = mech_b.update_species(self.lacI, self.assembly.dna, n_mer_species=self.lacI, component = self, part_id = self.assembly.dna.name+'_2x_'+self.lacI.name, cooperativity=2)
species += species_b
#print('plac stuff \n',species_b)
for s in species_b:
if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.assembly.dna in s.species:
self.plac_2_lacI = s
if self.leak is not False:
species += mech_tx.update_species(dna = self.plac_2_lacI, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
# IPTG can bind to lacI to sequester it
# give citation for 1 being enough
"""A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
2lacI <--> 2xlacI
2xlacI + IPTG <--> 2xlacI:IPTG
"""
species_b = mech_b.update_species(self.lacI, self.IPTG, component = self, part_id = self.IPTG.name+'_2x_'+self.lacI.name, cooperativity=2)
#print('iptg stuff \n', species_b)
species += species_b
for s in species_b:
if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.IPTG in s.species:
self.IPTG_2_lacI = s
# Unbound transcription
species += mech_tx.update_species(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
return species
def update_reactions(self):
print('r')
mech_tx = self.mechanisms['transcription']
mech_b = self.mechanisms['binding']
if self.IPTG_2_lacI == None or self.plac_2_lacI == None:
self.update_species()
if self.leak != False:
# Will need transcription, plac_2x_lacI, ktx, ku, kb
reactions += mech_tx.update_reactions(dna = self.plac_2_lacI, component = self, part_id = self.plac_2_lacI.name, \
transcript = self.transcript, protein = self.assembly.protein)
reactions = []
# Repression
# Will need two_step_cooperative_binding, plac_2x_lacI, kb1 ku1, kb2, ku2
reactions += mech_b.update_reactions(self.lacI, self.assembly.dna, component = self, \
part_id = self.plac_2_lacI.name, cooperativity = 2)
# LacI binding to IPTG
# Will need two_step_cooperative_binding, IPTG_2x_lacI, kb1 ku1, kb2, ku2
# .name gets rid of the complex
# .type will give you complex, species, etc
# .attributes gives you properties like phosphorylation state
reactions += mech_b.update_reactions(self.lacI, self.IPTG, component = self, part_id = self.IPTG_2_lacI.name, cooperativity=2)
'''
Any bit of IPTG destroys the bound structure
IPTG + plac_2xlacI --> IPTG:2xlacI + plac
'''
if isinstance(mech_b, Two_Step_Cooperative_Binding):
kb = self.get_parameter("kb2", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
else:
kb = self.get_parameter("kb", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
print(self.IPTG_2_lacI.name, kb)
reactions.append(Reaction([self.IPTG, self.plac_2_lacI], [self.IPTG_2_lacI, self.assembly.dna], k=kb))
# Unbound transcription plac, kb, ku, ktx
reactions += mech_tx.update_reactions(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
return reactions
Sort out parameter printer¶
In [43]:
def param_printer2(x0, parameters):
inits = pd.DataFrame(list(x0.items()),columns = ['param','value'])
df = pd.DataFrame(list(parameters.items()),columns = ['param','value'])
col_part_id = []
mech = []
for val in df['param']:
if isinstance(val, tuple):
col_part_id.append(val[0])
mech.append(val[1])
else:
col_part_id.append(' ')
mech.append(val)
df['param'] = mech
df['part'] = col_part_id
df.reindex(['part','param','value'],axis=1)
df[''] = ['|']*len(col_part_id)
return pd.concat([df, inits], axis=1).fillna("")
AND simulation¶
In [44]:
parameters = {
('pttr_phosphate_P_protein_ttrR', 'cooperativity'): 2.0,
('pttr_phosphate_P_protein_ttrR', 'ktx'): 0.17025, #this is the transcription rate
('pttr_phosphate_P_protein_ttrR', 'ku'): 11.01321586, #this is the polymerase unbinding rate
('transcription_mm', 'ktx'): .01, #These are the parameters for transcription leak
('transcription_mm','kb'):100, #default binding rate!
#('ttrS-P',"kdephos"):1000, #auto dephosphorylation
('ligand_tt_protein_ttrS-P','kdephos'):1000, #dephosphorylation rate of ttrS
#('ttrS-P',"kphos"):10, #ttrs phosphorylates itself automatically. This is a proxy for sensing
('ligand_tt_protein_ttrS-P','kphos'):10, #equivalent to kphos, from before
#('ttrS-ttrR','ktransfer'):5, #rate of transfer from ttrS to ttrR
('ligand_tt_protein_ttrS-ttrR', 'ktransfer'):5, #this is the new transfer rate
('ttrR-P',"kdephos"):1000, #auto dephosphorylation
"kb":100, "ku":10, "ktx":.05, "ktl":.2, "kdeg":3, 'cooperativity':2, #some default parameters
('phrpL_hrpR', 'cooperativity'): 1.0, #default is 2, but these proteins bind at cooperativity = 1
('phrpL_hrpS', 'cooperativity'): 1.0,
('phrpL_hrpR_hrpS_RNAP', 'ktx'): 0.17025, #this is the transcription rate
('phrpL_hrpR_hrpS_RNAP', 'ku'): 11.01321586, #this is the polymerase unbinding rate
('translation_mm', 'B0030', 'ku'): 888, #Unbinding
('translation_mm', 'B0030', 'ktl'): .15, #Translation Rate
('transcription_mm', 'ktx'): .01, #These are the parameters for transcription leak
('transcription_mm', 'ku'): 100,
('phrpL_hrpR','ku'):100,
('phrpL_hrpS','ku'):100,
('phrpL_hrpS_RNAP', 'ku'): 1000.,#we'll say that leak has 100x the unbinding rate
('phrpL_hrpR_RNAP', 'ku'): 1000,
('phrpL','ku'):100000,
('binding', 'pLac_lacI','ktx'):0.001,
('plac_lacI','kb'):10,
('plac_lacI','ku'):0.001,
('binding', 'plac_hrpS_2x_lacI','kb'):8888,
("simple_transcription", 'plac','ktx'):0.1,
#Make a lot of repressor
('translation_mm', 'rbsL', 'ku'): 0.0001, #Unbinding
('translation_mm', 'rbsL', 'ktl'): 2, #Translation Rate
('plac_hrpS','kb'):100,
('plac_hrpS','ku'):0.0001
}
#P_rbsS_ttrS
constit_ttrS = DNAassembly(name = "constit_ttrS", promoter = "P", rbs = "rbsS", protein = "ttrS")
#P_rbsR_ttrR
constit_ttrR = DNAassembly(name = "constit_ttrR", promoter = "P", rbs = "rbsR", protein = "ttrR")
#P_arl_lacI
constit_lacI = DNAassembly(name = "constit_lacI", promoter = "P", rbs = "rbsL", protein = "lacI")
#P_arl_lacI
constit_mScarlet = DNAassembly(name = "constit_mScarlet", promoter = "P", rbs = "B0034", protein = "mScarlet")
# plac_hrpS
IPTG = Species("IPTG", material_type='inducer')
#IPTG = Protein("IPTG")
lac_p = LacIPromoter2(name="pLac", IPTG = "IPTG", lacI = "lacI", leak=False)
plac_hrpS = DNAassembly(name = "plac_hrpS", promoter=lac_p, rbs="rbsH",protein="hrpS")
# pttr_hrpR
tt = Species("tt",material_type="ligand")
ttrSbound = ChemicalComplex([Protein("ttrS"),tt])
ttrS3 = PhosphoTransferase(ttrSbound.species,targets=["ttrR"])
ttrR = PhosphoProtein("ttrR")
# Need to change from UTR1
pttr_hrpR = DNAassembly(name = "pttr_hrpR",promoter=RegulatedPromoter("pttr",[ttrR.get_phosphorylated_species()],leak=False),rbs="UTR1",protein="hrpR")
hrpL_gfp = DNAassembly("hrpL_gfp",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=False),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli", components = [pttr_hrpR,ttrSbound,ttrS3,ttrR,constit_ttrS,constit_ttrR,
constit_lacI, plac_hrpS, hrpL_gfp, constit_mScarlet],parameters=parameters)
CRN_extract_1 = extract_1_TXTL.compile_crn()
In [45]:
param_printer2(x0, parameters)
Out[45]:
In [46]:
x0 = {"dna_pttr_hrpR":1,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
"dna_constit_lacI":1,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":1,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":15,
"protein_RNAase":30,
"protein_Ribo":120,
"protein_lacI": 5}
In [47]:
IPTG_um_list = np.logspace(-3, 2, 6).round(decimals=6)
tt_um_list = np.logspace(-3, 2, 6).round(decimals=6)
fig, heat_map = simulate_heatmap('ligand_tt',tt_um_list,'inducer_IPTG',IPTG_um_list,x0=x0,normalize=True,title='Compiled Full AND gate \n lacI Spike', CRN=CRN_extract_1)
heat_map.set(xlabel='IPTG (uM)', ylabel='Tetrathionate (uM)')
heat_map.invert_yaxis()
plt.yticks(fontsize=18, fontstyle='normal', rotation = 0)
fig.savefig("./figures/fig6/fig6d.png")
plt.show()
100) {
clearInterval(timer);
console.log("Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing");
}
}
}, 10, root)
}
})(window);
In [36]:
all_sims['Time (hr)'] = timepoints
all_sims.to_csv('./data/all_sims.csv', index = False)
HRP AND Circuit¶
This $\sigma 54$-dependent split activator requires two proteins, hrpR and hrpS, to function. For now, let's just spike in the two proteins and see how our and gate responds.
In [37]:
txtl = CRNLab("GFP")
txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = './parameters/BasicExtract.tsv')
#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
Re1.head()
Out[37]:
Vary hrpR and hrpS proteins¶
Now, we will make sure the AND gate (combinatorial promoter) is working by varying protein concentrations.
Get plotting in order¶
Let's create a plotting function. We will plot a heatmap of protein GFP at the last slice of the simulation. Because we have no degradation of protein, the last timepoint will be the max. We will be able to pass in normalize, which divides all values by the maximum of all simulated values, or max_to_1, which will set the scale bar to be 0 to 1. This basically causes plots that might've been
In [38]:
def simulate_heatmap(inducer_1_name,
inducer_1_values,
inducer_2_name,
inducer_2_values,
title,
x0=x0,
normalize=True,
max_to_1 = False,
CRN=CRN_extract_1,
timepoints=timepoints):
GFP_max = pd.DataFrame()
for conc_1 in inducer_1_values:
x0[inducer_1_name] = conc_1
for conc_2 in inducer_2_values:
x0[inducer_2_name] = conc_2 #Change my initial condition dictionary
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
GFP_max = GFP_max.append({inducer_1_name:conc_1,
inducer_2_name:conc_2,
'GFP_max': Re1["protein_GFP"].values[-1]}, ignore_index=True)
data = pd.pivot_table(data = GFP_max, index = inducer_1_name,
columns = inducer_2_name,
values = 'GFP_max')
fig, ax = plt.subplots(figsize=(6, 4.5))
if normalize:
data /= np.max(GFP_max['GFP_max'])
if max_to_1:
heat_map = sns.heatmap(data,
cmap = 'viridis',
linewidths=.1, vmin = 0, vmax = 1,
cbar_kws={'label': 'Normalized GFP'})
else:
heat_map = sns.heatmap(data,
cmap = 'viridis',
linewidths=.2,
cbar_kws={'label': 'Normalized GFP'})
sns.set(font_scale=1.4)
# plt.tick_params(axis='both', labelsize=20)
plt.xticks(fontsize=18, fontstyle='normal', rotation = 90)
plt.yticks(fontsize=18, fontstyle='normal')
plt.xlabel(inducer_1_name, fontsize=18, fontweight='bold')
plt.ylabel(inducer_2_name, fontsize=18, fontweight='bold')
heat_map.set_title(title)
fig.tight_layout()
return fig, heat_map
And Gate with leak in hrpR¶
Let's say hrpR mutated and now it can recruit RNAP better than it did beter. We will change phrpL_hrpR_RNAP ku = 5 instead of 500.
In [39]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak5_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 1x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_1.png", bbox_inches = "tight")
plt.show()
In [40]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak50_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 10x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_2.png", bbox_inches = "tight")
plt.show()
In [41]:
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 100x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_3.png", bbox_inches = "tight")
plt.show()
Full Compiled AND gate with expression of hrpR, S¶
LacI promoter class¶
In [42]:
class LacIPromoter2(Promoter):
def __init__(self, name = "pLac", IPTG = "IPTG",
lacI = "lacI", leak = True, assembly = None,
transcript = None, length = 0, mechanisms = {},
parameters = {}, **keywords):
Promoter.__init__(self, name=name, transcript=transcript, assembly=assembly, length=length, parameters=parameters, **keywords)
#RegulatedPromoter.__init__(self = self, name = name, regulators = regulators)
# Transcription is already included in promoter class
self.default_mechanisms = {"binding": One_Step_Cooperative_Binding()}
self.leak = leak
self.IPTG = self.set_species(IPTG, material_type = "inducer")
self.lacI = self.set_species(lacI, material_type = "protein")
self.plac_2_lacI = None
self.IPTG_2_lacI = None
def update_species(self):
mech_tx = self.mechanisms['transcription']
mech_b = self.mechanisms['binding']
species = []
self.complexes = []
"""A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
2lacI <--> 2xlacI
2xlacI + pLac <--> 2xlacI:pLac
"""
species_b = mech_b.update_species(self.lacI, self.assembly.dna, n_mer_species=self.lacI, component = self, part_id = self.assembly.dna.name+'_2x_'+self.lacI.name, cooperativity=2)
species += species_b
#print('plac stuff \n',species_b)
for s in species_b:
if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.assembly.dna in s.species:
self.plac_2_lacI = s
if self.leak is not False:
species += mech_tx.update_species(dna = self.plac_2_lacI, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
# IPTG can bind to lacI to sequester it
# give citation for 1 being enough
"""A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
2lacI <--> 2xlacI
2xlacI + IPTG <--> 2xlacI:IPTG
"""
species_b = mech_b.update_species(self.lacI, self.IPTG, component = self, part_id = self.IPTG.name+'_2x_'+self.lacI.name, cooperativity=2)
#print('iptg stuff \n', species_b)
species += species_b
for s in species_b:
if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.IPTG in s.species:
self.IPTG_2_lacI = s
# Unbound transcription
species += mech_tx.update_species(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
return species
def update_reactions(self):
print('r')
mech_tx = self.mechanisms['transcription']
mech_b = self.mechanisms['binding']
if self.IPTG_2_lacI == None or self.plac_2_lacI == None:
self.update_species()
if self.leak != False:
# Will need transcription, plac_2x_lacI, ktx, ku, kb
reactions += mech_tx.update_reactions(dna = self.plac_2_lacI, component = self, part_id = self.plac_2_lacI.name, \
transcript = self.transcript, protein = self.assembly.protein)
reactions = []
# Repression
# Will need two_step_cooperative_binding, plac_2x_lacI, kb1 ku1, kb2, ku2
reactions += mech_b.update_reactions(self.lacI, self.assembly.dna, component = self, \
part_id = self.plac_2_lacI.name, cooperativity = 2)
# LacI binding to IPTG
# Will need two_step_cooperative_binding, IPTG_2x_lacI, kb1 ku1, kb2, ku2
# .name gets rid of the complex
# .type will give you complex, species, etc
# .attributes gives you properties like phosphorylation state
reactions += mech_b.update_reactions(self.lacI, self.IPTG, component = self, part_id = self.IPTG_2_lacI.name, cooperativity=2)
'''
Any bit of IPTG destroys the bound structure
IPTG + plac_2xlacI --> IPTG:2xlacI + plac
'''
if isinstance(mech_b, Two_Step_Cooperative_Binding):
kb = self.get_parameter("kb2", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
else:
kb = self.get_parameter("kb", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
print(self.IPTG_2_lacI.name, kb)
reactions.append(Reaction([self.IPTG, self.plac_2_lacI], [self.IPTG_2_lacI, self.assembly.dna], k=kb))
# Unbound transcription plac, kb, ku, ktx
reactions += mech_tx.update_reactions(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
return reactions
Sort out parameter printer¶
In [43]:
def param_printer2(x0, parameters):
inits = pd.DataFrame(list(x0.items()),columns = ['param','value'])
df = pd.DataFrame(list(parameters.items()),columns = ['param','value'])
col_part_id = []
mech = []
for val in df['param']:
if isinstance(val, tuple):
col_part_id.append(val[0])
mech.append(val[1])
else:
col_part_id.append(' ')
mech.append(val)
df['param'] = mech
df['part'] = col_part_id
df.reindex(['part','param','value'],axis=1)
df[''] = ['|']*len(col_part_id)
return pd.concat([df, inits], axis=1).fillna("")
AND simulation¶
In [44]:
parameters = {
('pttr_phosphate_P_protein_ttrR', 'cooperativity'): 2.0,
('pttr_phosphate_P_protein_ttrR', 'ktx'): 0.17025, #this is the transcription rate
('pttr_phosphate_P_protein_ttrR', 'ku'): 11.01321586, #this is the polymerase unbinding rate
('transcription_mm', 'ktx'): .01, #These are the parameters for transcription leak
('transcription_mm','kb'):100, #default binding rate!
#('ttrS-P',"kdephos"):1000, #auto dephosphorylation
('ligand_tt_protein_ttrS-P','kdephos'):1000, #dephosphorylation rate of ttrS
#('ttrS-P',"kphos"):10, #ttrs phosphorylates itself automatically. This is a proxy for sensing
('ligand_tt_protein_ttrS-P','kphos'):10, #equivalent to kphos, from before
#('ttrS-ttrR','ktransfer'):5, #rate of transfer from ttrS to ttrR
('ligand_tt_protein_ttrS-ttrR', 'ktransfer'):5, #this is the new transfer rate
('ttrR-P',"kdephos"):1000, #auto dephosphorylation
"kb":100, "ku":10, "ktx":.05, "ktl":.2, "kdeg":3, 'cooperativity':2, #some default parameters
('phrpL_hrpR', 'cooperativity'): 1.0, #default is 2, but these proteins bind at cooperativity = 1
('phrpL_hrpS', 'cooperativity'): 1.0,
('phrpL_hrpR_hrpS_RNAP', 'ktx'): 0.17025, #this is the transcription rate
('phrpL_hrpR_hrpS_RNAP', 'ku'): 11.01321586, #this is the polymerase unbinding rate
('translation_mm', 'B0030', 'ku'): 888, #Unbinding
('translation_mm', 'B0030', 'ktl'): .15, #Translation Rate
('transcription_mm', 'ktx'): .01, #These are the parameters for transcription leak
('transcription_mm', 'ku'): 100,
('phrpL_hrpR','ku'):100,
('phrpL_hrpS','ku'):100,
('phrpL_hrpS_RNAP', 'ku'): 1000.,#we'll say that leak has 100x the unbinding rate
('phrpL_hrpR_RNAP', 'ku'): 1000,
('phrpL','ku'):100000,
('binding', 'pLac_lacI','ktx'):0.001,
('plac_lacI','kb'):10,
('plac_lacI','ku'):0.001,
('binding', 'plac_hrpS_2x_lacI','kb'):8888,
("simple_transcription", 'plac','ktx'):0.1,
#Make a lot of repressor
('translation_mm', 'rbsL', 'ku'): 0.0001, #Unbinding
('translation_mm', 'rbsL', 'ktl'): 2, #Translation Rate
('plac_hrpS','kb'):100,
('plac_hrpS','ku'):0.0001
}
#P_rbsS_ttrS
constit_ttrS = DNAassembly(name = "constit_ttrS", promoter = "P", rbs = "rbsS", protein = "ttrS")
#P_rbsR_ttrR
constit_ttrR = DNAassembly(name = "constit_ttrR", promoter = "P", rbs = "rbsR", protein = "ttrR")
#P_arl_lacI
constit_lacI = DNAassembly(name = "constit_lacI", promoter = "P", rbs = "rbsL", protein = "lacI")
#P_arl_lacI
constit_mScarlet = DNAassembly(name = "constit_mScarlet", promoter = "P", rbs = "B0034", protein = "mScarlet")
# plac_hrpS
IPTG = Species("IPTG", material_type='inducer')
#IPTG = Protein("IPTG")
lac_p = LacIPromoter2(name="pLac", IPTG = "IPTG", lacI = "lacI", leak=False)
plac_hrpS = DNAassembly(name = "plac_hrpS", promoter=lac_p, rbs="rbsH",protein="hrpS")
# pttr_hrpR
tt = Species("tt",material_type="ligand")
ttrSbound = ChemicalComplex([Protein("ttrS"),tt])
ttrS3 = PhosphoTransferase(ttrSbound.species,targets=["ttrR"])
ttrR = PhosphoProtein("ttrR")
# Need to change from UTR1
pttr_hrpR = DNAassembly(name = "pttr_hrpR",promoter=RegulatedPromoter("pttr",[ttrR.get_phosphorylated_species()],leak=False),rbs="UTR1",protein="hrpR")
hrpL_gfp = DNAassembly("hrpL_gfp",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=False),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli", components = [pttr_hrpR,ttrSbound,ttrS3,ttrR,constit_ttrS,constit_ttrR,
constit_lacI, plac_hrpS, hrpL_gfp, constit_mScarlet],parameters=parameters)
CRN_extract_1 = extract_1_TXTL.compile_crn()
In [45]:
param_printer2(x0, parameters)
Out[45]:
In [46]:
x0 = {"dna_pttr_hrpR":1,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
"dna_constit_lacI":1,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":1,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":15,
"protein_RNAase":30,
"protein_Ribo":120,
"protein_lacI": 5}
In [47]:
IPTG_um_list = np.logspace(-3, 2, 6).round(decimals=6)
tt_um_list = np.logspace(-3, 2, 6).round(decimals=6)
fig, heat_map = simulate_heatmap('ligand_tt',tt_um_list,'inducer_IPTG',IPTG_um_list,x0=x0,normalize=True,title='Compiled Full AND gate \n lacI Spike', CRN=CRN_extract_1)
heat_map.set(xlabel='IPTG (uM)', ylabel='Tetrathionate (uM)')
heat_map.invert_yaxis()
plt.yticks(fontsize=18, fontstyle='normal', rotation = 0)
fig.savefig("./figures/fig6/fig6d.png")
plt.show()
In [48]:
x0 = {"dna_pttr_hrpR":1,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
"dna_constit_lacI":1,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":1,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":15,
"protein_RNAase":30,
"protein_Ribo":120,
"protein_lacI": 0}
In [49]:
IPTG_um_list = np.logspace(-3, 2, 6).round(decimals=6)
tt_um_list = np.logspace(-3, 2, 6).round(decimals=6)
fig, heat_map = simulate_heatmap('ligand_tt',tt_um_list,'inducer_IPTG',IPTG_um_list,x0=x0,normalize=True,title='Compiled Full AND gate \n No lacI Spike', CRN=CRN_extract_1)
heat_map.set(xlabel='IPTG (uM)', ylabel='Tetrathionate (uM)')
heat_map.invert_yaxis()
plt.yticks(fontsize=18, fontstyle='normal', rotation = 0)
fig.savefig("./figures/fig6/fig6b.png")
plt.show()
Early Circuit Dynamics in no lacI Spike Model¶
In [50]:
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_pttr_hrpR":5,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
# Cole origin
"dna_constit_lacI":40,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":1,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":10,
"protein_Ribo":90.,
"ligand_tt": 1,
"inducer_IPTG": 0.,
"protein_lacI": 0,}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
results_melted = pd.melt(Re1, var_name='species', id_vars=['time'])
results_melted.head()
Out[50]:
In [51]:
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
In [52]:
p = bokeh.plotting.figure(title = 'Early Circuit Dynamics, pLac', x_axis_label = 'Time (min)', y_axis_label = 'Simulated Concentration (uM)', height = 2000, width = 4300)
p.line(timepoints, Re1['complex_dna_plac_hrpS_protein_RNAP'], legend_label = 'Active pLac', color = 'Green', line_width = 14)
p.line(timepoints, Re1['rna_constit_lacI'], legend_label = 'LacI RNA', color = 'orange', line_width = 14)
p.line(timepoints, Re1['protein_hrpS'], legend_label = 'hrpS protein', color='lightsteelblue', line_width = 14)
p.line(timepoints, Re1['complex_dna_hrpL_gfp_protein_hrpR_protein_hrpS_protein_RNAP'], legend_label = 'hrpS, hrpR both bound to pHRPL', color='darkorchid', line_width = 14)
p.line(timepoints, Re1['rna_hrpL_gfp'], legend_label = 'GFP RNA', color='navy', line_width = 14)
p.add_layout(p.legend[0], 'right')
p.legend.title = 'Species'
p.legend.glyph_width = 100
p.legend.padding = 100
p.legend.title_text_font_size = '80pt'
p.legend.label_standoff = 50
p.legend.glyph_height = 100
p.xaxis.major_tick_line_width = 10
p.yaxis.major_tick_line_width = 10
p.xaxis.minor_tick_line_width = 10
p.yaxis.minor_tick_line_width = 10
p.axis.minor_tick_out = 20
p.axis.minor_tick_out = 20
p.axis.major_tick_out = 40
p.axis.major_tick_out = 40
p.yaxis.axis_label_standoff = 50
p.xaxis.axis_label_standoff = 50
p.xaxis.axis_line_width = 8
p.yaxis.axis_line_width = 8
p.yaxis.major_label_standoff = 30
p.xaxis.major_label_standoff = 30
p.ygrid.grid_line_width = 10
p.xgrid.grid_line_width = 10
p.axis.axis_label_text_font_size = "80pt"
p.legend.label_text_font_size = "80pt"
p.title.text_font_size = '100pt'
p.yaxis.major_label_text_font_size = "80pt"
p.xaxis.major_label_text_font_size = "80pt"
#p.y_range = Range1d(0, 0.9)
bokeh.io.show(p)
Computing Environment¶
In [53]:
%load_ext watermark
%watermark -v -p numpy,pandas,scipy,bioscrape,biocrnpyler,bokeh,jupyterlab,holoviews